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Persistence probabilities of the interface height in (1+1)- and (2+1) -dimensional atomistic, solid- 
on-solid, stochastic models of surface growth are studied using kinetic Monte Carlo simulations, with 
emphasis on models that belong to the molecular beam epitaxy (MBE) universality class. Both the 
initial transient and the long-time steady-state regimes are investigated. We show that for growth 
models in the MBE universality class, the nonlinearity of the underlying dynamical equation is 
clearly reflected in the difference between the measured values of the positive and negative persistence 
exponents in both transient and steady-state regimes. For the MBE universality class, the positive 
and negative persistence exponents in the steady-state are found to be 9% = 0.66 ± 0.02 and 9l = 
0.78+0.02, respectively, in (1+1) dimensions, and Of = 0.76+0.02 and 6^ = 0.85+0.02, respectively, 
in (2+1) dimensions. The noise reduction technique is applied on some of the (l+l)-dimensional 
models in order to obtain accurate values of the persistence exponents. We show analytically that 
a relation between the steady-state persistence exponent and the dynamic growth exponent, found 
earlier to be valid for linear models, should be satisfied by the smaller of the two steady-state 
persistence exponents in the nonlinear models. Our numerical results for the persistence exponents 
are consistent with this prediction. We also find that the steady-state persistence exponents can be 
obtained from simulations over times that are much shorter than that required for the interface to 
reach the steady state. The dependence of the persistence probability on the system size and the 
sampling time is shown to be described by a simple scaling form. 

PACS numbers: 81.10.Aj, 05.20.-y, 05.40.-a, 68.35.Ja 



I. INTRODUCTION 



Nonequilibrium surface growth and interface dynamics represent an area of research that has received much atten- 
tion in the last two decades 01 . A larg e number of discrete atomistic growth models [E IE EL IE IE IE 13 an d stochastic 
growth equations |ij Hoi Ull O EE El have been found [3 to exhibit generic scale invariance characterized by power 
law behavior of several quantities of interest, such as the interface width as a function of time (measured in units of 
deposited layers) and space- and time-dependent correlation functions of the interface height. Much effort has been 
devoted to the classification of growth models and equations into different universality classes characterized by the 
values of the exponents that describe the dynamic scaling behavior implied by these power laws. A variety of experi- 
mental studies [l5lll6j| have confirmed the occurrence of dynamic scaling in nonequilibrium epitaxial growth. Among 
the various experimental methods of surface growth, molecular beam epitaxy (MBE) is especially important because 
it plays a crucial role in the fabrication of smooth semiconductor films required in technological applications. Under 
usual MBE growth conditions, desorption from the film surface is negligible and the formation of bulk vacancies and 
overhangs is strongly suppressed, ft is generally believed that nonequilibrium surface growth under these conditions is 
well-described by a conserved nonlinear Langevin-type equation 0, 0, 0] and related atomistic models |E IE IE Q 
that form the so-called "MBE universality class" . 

Surface growth is an example of a general class of problems involving the dynamics of non-Markovian, spatially 
extended, stochastic systems. In recent years, t he conc ept of p ersistence [l7| has proven to be very useful in analyzing 
the dynamical behavior of such systems EE EE HE EB 113 uE Q • Loosely speaking, a stochastic variable is persistent 
if it has a tendency to maintain its initial characteristics over a long period of time. The persistence probability P(t) 
is typically defined as the probability that a characteristic feature (e.g. the sign) of a stochastic variable does not 
change at all over a certain period of time t. Although the mathematical concept of persistence was introduced a 
long time ago in the context of the "zero-crossing problem" in Gaussian stationary processes |25j |. it is only very 
recently that this concept has received attention in describing the statistics of first passage events in a variety of 
spatially extended nonequilibrium systems. Examples of such applications of the concept of persistence range from 
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the fundamental classical diffusion equation [18( to the zero temperature Glauber dynamics of the ferromagnetic Ising 
and g-state Potts models [l|l |2(J, W& l2ij and phase ordering kinetics [23 ■ Recently, a generalization of the persistence 
concept (probability of persistent l arge deviations) has been introduced (2|| . A closely related idea, that of sign-time 
distribution, was developed in Ref. |2jj • An increasing number of experimental results are also available for persistence 
in systems such as coalescence of droplets [28| , coarsening of two-dimensional soap froth [2!| , twisted nematic liquid 
crystal [23, and nuclear spin distribution in laser polarized Xe 129 gas [3lj . 

Recent work of Krug and collaborators [2^, 0] has extended the persistence concept to the first-passage statistics 
of fluctuating interfaces. Persistence in the dynamics of fluctuating interfaces is of crucial importance in ultra-small 
scale solid-state devices. As the technology advances into the nanometric regime, questions such as how long a 
particular perturbation that appears in an evolving interface persists in time and what is the average time required 
for a structure to first fluctuate into an unstable configuration become im por tant. The persistence probability can 
provide quantitative predictions on such questions. Recent experiments |.32i 18.1 1.3 1| have demonstrated the usefulness 
of the concept of persistence in the characterization of the equilibrium fluctuations of ste ps o n a vicinal surface. 
Analysis of experimental data on step fluctuations on Al/Si(lll) [2^, [3] and Ag(lll) [33l l34| surfaces has shown 
that the long-time behavior of the persistence probability and the probability of persistent large deviations in these 
systems agrees quantitatively with the corresponding theoretical predictions. These results show that the persistence 
probability and related quantities are particularly relevant for describing and understanding the long time dynamics 
of interface fluctuations. 

In the context of surface growth and fluctuations, the persistence probability P(t ,t + t) may be defined as the 
probability that starting from an initial time to, the interfacial height h(r,t') at spatial position r does not return 
to its original value at any point in the time interval between t and t + t. This probability is clearly the sum of 
the probabilities of the height h(r,t') always remaining above (the positive persistence probability P + ) and always 
remaining below (the negative persistence probability P_) its specific initial value h(r,to) for all to < t' < to + t. This 
concept quantifies the tendency of a stochastic field (in our case the interface height) to persistently conserve a specific 
feature (the sign of the interfacial height fluctuations). The persistence probability P(to,to + t) would, in general, 
depend on both to and t. In the early stage of the growth process starting from a flat interface (transient regime), the 
interface gradually develops dynamical roughness [l| due to the effect of fluctuations in the beam intensity. In this 
regime, the choice of the initial time to is clearly important: it determines the degree of roughness of the configuration 
from which the interface evolves. At long times, the growing interface enters into a new evolution stage, called the 
steady-state regime, characterized by fully developed roughness that does not increase further in time. In this regime, 
the choice of to is expected to be unimportant. 

The work of Krug et al. [23| shows that for a class of linear Langevin-type equations for surface growth and atomistic 
models belonging in the same dynamical universality class as these equations, the persistence probability decays as 
a power law in time for long times in both transient and steady-state regimes. These power laws define the positive 
and negative persistence exponents, 9± and #±, for positive and negative persistence in the transient and steady-state 
regimes, respectively. The h — > —h symmetry of the linear growth equations implies that 9+ = (PL and 9+ = 9^_ in 
these systems. In Ref. [23, it was pointed out that the persistence exponent in the steady state of these linear models 
is related to the dynamic scaling exponent (3, that describes the growth of the interface width W as a function of time 
t in the transient regime (W cx f 9 ), through the relation 9+ = 9^ = 1 — (3. The validity of this relation was confirmed 
by numerical simulations. Since the exponent (3 is the same for all models in the same dynamical universality class, 
this result implies that the persistence exponent in the steady-state regime of these linear models is also universal. 
Numerical results for the persistence exponent in the transient regime, for which no analytic predictions are available, 
also indicate a similar universality. Kallabis and Krug |24j carried out a similar calculation for (l+l)-dimensional 
Kardar-Parisi-Zhang (KPZ) [ll| interfaces. They found that the nonlinearity in the KPZ equation that breaks the 
h — > —h symmetry is reflected in different values of the positive and negative persistence exponents, and in the 
transient regime. The values of the steady-state persistence exponents 9+ and 9^_ were found to be equal to each other, 
and equal to 1 — (3 within the accuracy of the numerical results. This is expected because the h — > —h symmetry is 
dynamically restored in the steady state of the (l+l)-dimensional KPZ equation. This is, however, a specific feature 
of the (l+l)-dimcnsional KPZ model, which for nongeneric reasons, turns out to be up-down symmetric in the steady 
state. Nonlinear surface growth models (e.g. the higher dimensional KPZ model, the nonlinear MBE growth model) 
are generically expected to have different values of 9± in both transient and steady-state regimes. 

In this paper, we present the results of a detailed numerical study of the persistence behavior of several atomistic, 
solid-on-solid (SOS) models of surface growth in (1+1) and (2+1) dimensions. While we concentrate on models 
in the MBE universality class, results for a few other models, some of which have been studied in Refs. [2^ and 
I24I are also presented for completeness. The highly non-trivial nature of the persistence probability, in spite of a 
deceptive simplicity of the defining concept, arises from the complex temporal non-locality ("memory") inherent 
in its definition. In fact, there are very few stochastic problems where an analytical solution for the persistence 
probability has been achieved. These include the classical Brownian motion [3^, the random acceleration problem 



3 



[36( and the one dimensional Ising and g-state Potts models 20]. In general, the highly nonlocal nature of the 
temporal correlations in a non-Markovian stochastic process makes it extremely difficult to obtain exact results for 
the persistence probability even for seemingly simple stochastic processes. Even for the simple diffusion equation, 
the persistence exponent is known only numerically, or within an independent interval approximation |l8j or series 
expansion approach |37| . However, it is fairly straightforward in most cases to directly simulate the persistence 
probability to obtain its stationary power law behavior at large times, and thus to numerically obtain the approximate 
value of the persistence exponent. For this reason, we use stochastic (Monte Carlo) simulations of the atomistic growth 
models to study their temporal persistence behavior in the transient and steady-state regimes. These models are 
defined in terms of random deposition and specific cellular-automaton-type local diffusion or relaxation rules. Some 
of these models are of the "limited-mobility" type in the sense that the surface diffusion rules or local restrictions limit 
the characteristic length over which a deposited particle can diffuse to just one or a few lattice spacings. The models 
in the MBE universality class considered in our study are: the Das Sarma-Tamborenea model j3j, the Wolf- Villain 
model 0, |3£|, the Kim-Das Sarma model Q and its "controlled" version 0] and the restricted solid-on-solid 
(RSOS) model of Kim at al. Q. We also present results for the Family model Q that is known to belong to the 
Edwards- Wilkinson universality class and the restricted solid-on-solid (RSOS) model of Kim and Kosterlitz Q 
that is in the KPZ universality class. 

The main objective of our study is to examine the effects of the nonlinearity in the MBE growth equation on 
the persistence behavior. Unlike the (l+l)-dimcnsional KPZ equation, the nonlinearity in the MBE growth equation 
persists in the steady state in the sense that the height profile exhibits a clear asymmetry between the positive and 
negative directions (above and below the average height). Therefore, the positive and negative persistence exponents 
0+ and 9^_ are expected to have different values in these models. If this is the case, then the relation between the 
steady-state persistence exponent and the dynamic scaling exponent (3 found in linear models can not be valid for 
both Of and f?£, indicating that at least one of these exponents is a new, nontrivial one not related to the usual 
dynamic scaling exponents. The values of 6+ and 9^_ and their relation to j3, as well as the values of the transient 
persistence exponents 9+ and #T are the primary questions addressed in our study. We also investigate the universality 
of these exponents by measuring them for several models that are known to belong in the same universality class as 
far as their dynamic scaling behavior is concerned. To obtain accurate values of the exponents, the "noise reduction" 
technique [3!j is employed in some of the simulations of (l+l)-dimensional models. We also address some questions 
related to the methodology of calculating persistence exponents from simulations. Since the value of the dynamical 
exponent z is relatively large for models in the MBE universality class, the time required for reaching the steady state 
grows quickly as the sample size L is increased (t sat oc L z ). As a result, it is difficult to reach the steady state in 
simulations for large L. It is, therefore, useful to find out whether the value of the steady-state persistence exponents 
can be extracted from calculations of P(to,to + *) with to t sat oc L z . Another issue in this context involves the 
effects of the finitcness of the sample size L and the sampling time St (the time interval between two successive 
measurements of the height profile) on the calculated persistence probability. An understanding of these effects is 
needed for extracting reliable values of the persistence exponents from simulations that always involve finite values of 
L and St. Understanding the effects of L and St on the persistence analysis is not only important for our simulations, 
but is also important in the experimental measurements of persistence which invariably involve finite system size and 
sampling time. 

The main results of our study are as follows. We find that the positive and negative steady-state persistence expo- 
nents for growth models in the MBE universality class are indeed different from each other, reflecting the asymmetry 
of the interface arising from the presence of nonlinearities in the underlying growth equation. Our results for these 
exponent values are: Of = 0.66 ± 0.02 and #£ = 0.78 ± 0.02, respectively, in (1+1) dimensions; 9f = 0.76 ± 0.02 and 
= 0.85 ± 0.02 in (2+1) dimensions. The values of the positive and negative persistence exponents for different 
models are clearly correlated with the asymmetry of the "above" and " below" (defined relative to the mean interface 
height) portions of the interface. We show analytically that the smaller one of the two steady-state persistence expo- 
nents should be equal to (1 — (3). Thus, the relation = 1-/3 derived in Ref. [23j for linear surface growth models 
is expected to be satisfied by Of for the nonlinear models considered here. Our numerical results are consistent with 
this expectation: we find that the positive persistence exponent is indeed close to (1 — (3), while the negative one is 
significantly higher. Similar asymmetry is found for the persistence exponents in the transient regime with 9\ < 9^ 
in MBE growth. Within the uncertainties in the numerically determined values of the exponents, they are universal 
in the sense that different models in the same dynamic universality class yield very similar values for these exponents. 
For the models in the Edwards- Wilkinson and KPZ universality classes, we find results in agreement with those of 
earlier studies pLE^O] . 

Our simulations also reveal that a measurement of the steady-state persistence exponents is possible from simu- 
lations in which the initial time to is much smaller than the time (~ L z ) required for the interfacial roughness to 
saturate. A similar result was reported in Ref. |23J where it was found that the steady-state persistence exponent may 
be obtained from a calculation of P(to, to + 1) with t <C to <C L z . We find that the restriction t <C to is not necessary 
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for seeing a power law behavior of P±(to,to +t) - a power law with the steady-state exponents is found even if t 
is close to or somewhat larger than i . We exploit this finding in some of our persistence simulations for (2+1) - 
dimensional growth models which are more relevant to experiments. These results, however, also imply that it would 
be extremely difficult to measure the transient persistence exponents from real surface growth experiments. Finally, 
we show that the dependence of the steady-state persistence probability on the sample size L and the sampling time 
St is described by a simple scaling function of the variables t/L z and 5t/L z . This scaling description is similar to that 
found recently |40] for a different "persistence probability" , the survival probability, that measures the probability of 
the height not returning to its avera ge v alue (rather than the initial value) over a certain period of time. Although 
the "persistence" and the "survival" [4(j probability seem to be qualitatively similar in their definitions, the two are 
mathematically quite unrelated, and in fact, no exponent can be defined for the survival probability. In this paper 
we only discuss the persistence probability and the persistence exponent for surface growth processes. 

The rest of the paper is organized as follows: in Sec. Ill Al we briefly discuss the main universality classes and their 
corresponding dynamic equations and scaling exponents relevant for surface growth phenomena. Section fll Bl contains 
a short overview of the persistence probability concept from the interface fluctuations perspective. The discrete 
stochastic SOS growth models considered in our study are described in Sec. 1 1 1 1 1 In addition, we briefly describe 
in this section the noise reduction technique which is employed in some of our simulations. Section IIVI contains a 
detailed description of our main results: in Sec. IIV Al our (1+1) -dimensional simulation results for the transient 
and steady-state persistence exponents are presented, focusing mostly on models described by the nonlinear MBE 
dynamical equation. Section IIVBI contains the analytic derivation of a relation between the smaller steady-state 
persistence exponent and the dynamic growth exponent. In Sec. IIV CI we introduce an alternative approach for 
measuring the steady-state persistence exponents, using a relatively short "equilibration time" that is much shorter 
than the time required for reaching the true steady state. Section liV Dl contains the results of our (2+l)-dimensional 
persistence calculation for a selection of linear and nonlinear models. Simulation results that establish a scaling form 
of the dependence of the persistence probability on the sample size and the sampling time are presented in Sec. IIV El 
The final Sec. contains a summary of our main results and a few concluding remarks. 



The dynamic scaling behavior of stochastic growth equations may be classified into several universality classes. Each 
universality class is characterized by a set of scaling exponents [l| which depend on the dimensionality of the problem. 
These exponents are (a, (3, z), where a is the roughness exponent describing the dependence of the amplitude of height 
fluctuations in the steady-state regime (t ^> L z ) on the sample size L, (3 is the growth exponent that describes the 
initial power law growth of the interface width in the transient regime (1 « ( < Jj z ), and z is the dynamical exponent 
related to the system size dependence of the time at which the interface width reaches saturation. Note that z = % 
for all the models considered in this paper. To describe the interface evolution we use the single-valued function, 
h(r,t), which represents the height of the growing sample at position r and deposition time t. The interfacial height 
fluctuations are described by the root-mean-squared height deviation (or interface width) which is a function of the 
substrate size L and deposition time t: 



where h(t) is the average sample thickness. The width W(L,t) scales as W(L,t) oc t° for t <C L z and W(L,t) oc L a 
for t ^> L z L z being the equilibration time of the interface, when its stationary roughness is fully developed. 

Since it is convenient to write the evolution equations in terms of the deviation of the height from its spatial average 
value, h(r, t) — h(t), from now on we will denote by /i(r, t) the interface height fluctuation measured from the average 
height. Extensive studies of dynamic scaling in kinetic surface roughening (for an extended review see Ref. p"5jp have 
revealed the existence of (at least) four universality classes that are described, in the long wavelength limit, by the 
following continuum equations and sets of scaling exponents (a, (3, z), shown for the l+l(2+l)-dimensional cases, 
respectively: 

(1) The Edwards-Wilkinson (EW) second order linear equation: i, j, 2 (0 (log),0 (log), 2) 



II. STOCHASTIC GROWTH EQUATIONS AND PERSISTENCE PROBABILITIES 



A. Growth equations and dynamic scaling 



W{L,t) = ((h(v,t)-h{t)) 2 ) 1/2 , 



(1) 




(2) 



(3) 
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(3) The Mullins-Herring (MH) fourth order linear equation: §, |, 4 (1, |, 4) 

dh(r,t) 



dt 
and 

(4) The MBE fourth order nonlinear equation: ~ 1, — \, ~3(~|, —\, 

dh(r,t) 4 V72I/V7;. ^ + m2 



i/ 4 V 4 Mr,t)+»?(r,t), (4) 



= -^4V 4 ^(r,i) + A 22 V 2 |(V/i(r,t)| 2 +77(r,i), (5) 



where (i=2, 4) and Xj (j=2, 22) are constant. The quantity n(r,t) represents the noise term which accounts for 
the random fluctuations in the deposition rate. We assume that the noise has Gaussian distribution with zero mean 
and correlator 



<77(n, i 1 )7 7 (r 2 , t 2 )) - £W(n - r 2 )5(ti - h), (6) 

D being a constant related to the strength of the bare noise. Note that we do not include the (trivial) constant 
external deposition flux term in the continuum growth equations since that is easily eliminated by assuming that the 
height fluctuation h is always measured with respect to the average interface which is growing at a constant rate. 

The concepts of universality classes and scaling exponents have been widely used in the literature to analyze the 
kinetics of surface growth and fluctuations. Our study based on persistence probabilities is motivated by the possibility 
that the concept of persistence may provide an additional (and complementary) tool to analyze the surface growth 
kinetics. It addresses fundamental questions such as: is persistence an independent (and new) conceptual tool for 
studying surface fluctuations or essentially equivalent (or perhaps complementary) to dynamic scaling? and does 
persistence lead to the definition of new universality classes on the basis of the values of the persistence exponent? To 
answer these questions, we consider, for each of the four universality classes mentioned above (i.e. Eqs. (0-©), at 
least one growth model and investigate how the associated persistence exponents are related to the dynamic scaling 
exponents mentioned above. 



B. Transient and steady— state persistence probabilities 



Our goal is to calculate the positive and negative persistence probabilities (P±(to, to+t)) for a growing (fluctuating) 
interface in the transient and steady-state regimes. Here to is the initial time, and we are interested in evaluating the 
probability of the height at a fixed position remaining persistently above (P+) or below (P_ ) its initial value (i.e. its 
value at to by definition) during the time period between t and to + 1. If one considers the special case to — 0, when 
the interface is completely flat, then the quantity of interest is the probability that the interfacial height (measured 
from its spatial average) does not return to its initial zero value up to time t. This case is known as the transient (T) 
regime. For values of t that are small compared to the time scale for saturation of the interface width (t sat (L) oc L z ), 
the persistence probabilities in this regime are expected to exhibit a power law decay in time: 



P£(0,i)ocf±j , (7) 

where 8± are called the transient positive and negative persistence exponents. In the particular case of linear contin- 
uum growth equations, these exponents are equal because the symmetry under a change of sign of h(r, t) remains valid 
at all stages of the growth process. However, in the case of dynamics governed by nonlinear continuum equations, 
the lack of this "up-down" interfacial symmetry implies that P+ and P_ (and therefore, the exponents 0\ and QL) 
would, in general, be different from each other. No universal relationship between the transient positive and negative 
persistence exponents and the dynamic scaling exponents is known to exist for any one of the four universality classes 
mentioned above. 

On the other hand, if one considers to larger than t sa t{L), then the quantity of interest is the probability that the 
interfacial height at a fixed position does not return to its specific value at initial time to during the subsequent time 
interval between to and to + t. Instead of being flat, the interface morphology at time to has completely developed 
roughness, which produces persistence exponents that are different from the transient exponents defined earlier. This 
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case is known as the steady-state (S) regime, If i <C L z , one expects to obtain in this regime the steady-state 
persistence probability with a power law decay in time p3| 



where 9± are the steady-state positive and negative persistence exponents. It has been pointed out by Krug et al. 
[2^ | that for systems described by linear Langevin equation, the steady-state persistence exponents are related to the 
dynamic scaling exponent (3 in the following way: 



The exponent (3 is well known for linear Langevin equations for surface growth dynamics, and is given in d-dimensions 
by P = (1 — d/z)/2 for nonconserved white noise (Eq. ©), where z, the dynamical exponent, is here precisely equal 
to the power of the gradient operator entering the linear continuum dynamical growth equation [i.e. z = 2 in Eq. J2J; 
z = 4 in Eq. (0J]. The relation defined by Eq. @ holds true for the Langevin equations of Eqs. (J2J) and (@J, which are 
obviously linear, as well as for the special case of the (l+l)-dimensional KPZ equation of Eq. © 24], which, despite 
its nonlincarity, behaves as the linear EW equation in the steady state. Since the positive and negative exponents are 
expected to be different for general nonlinear Langevin equations, the relation of Eq. 10 can not be valid for both 
9+ and 9^ in systems described by such nonlinear equations. Therefore, at least one (or perhaps both) of these two 
persistence exponents must be non-trivial in the sense that it is not related to the usual dynamic scaling exponents. 
For this reason we pay particular attention to the MBE nonlinear equation and investigate whether its persistence 
exponents can be related to the dynamic scaling exponents. 



In this paper, we use different atomistic limited-mobility growth models for simulating surface growth processes. 
In these models, the substrate consists of a collection of lattice sites labeled by the index j (j = 1, 2, . . . , L d ) and the 
height variables h(xj) take integral values. The term "limited-mobility" is meant to imply that in these models, each 
adatom is characterized by a finite diffusion length which is taken to be one lattice spacing in most of the models we 
consider here. Thus, a deposited atom can explore only a few neighboring lattice sites according to a set of specific 
mobility rules before being incorporated into the growing film. The solid-on-solid constraint is imposed in all these 
models, so that defects such as overhangs and bulk vacancies are not allowed. In most of the models considered in this 
work, the possibility of desorption is neglected, thereby making the models "conserved" in the sense that all deposited 
atoms are incorporated in the film; the noise (given by Eq. ©) is of course nonconserved since the system is open to 
the deposition flux. 

The deposition process is described by a few simple rules in these models. An atomic beam drops atoms on 
the substrate in a random manner. Once a lattice site on the substrate is randomly chosen, the diffusion rules of 
the model are applied to the atom dropped at the chosen site to determine where it should be incorporated. The 
allocated site is then instantaneously filled by the adatom. We consider both (1+1)- and (2+1)- dimensional models 
(one or two spatial dimensions and one temporal dimension) defined on substrates of length L in units of the lattice 
spacing. The deposition rate is taken to be constant and equal to L d particles per unit time in our simulations of 
the Family (F), larger curvature (LC), Das Sarma-Tamborenea (DT), Wolf- Villain (WV) and controlled Kim-Das 
Sarma (CKD) models (see below). In these simulations, one complete layer is grown in each unit of time. In the 
RSOS Kim-Kosterlitz (KK) and Kim-Park-Kim (KPK) models described below, the diffusion rules are replaced by 
a set of local restrictions on nearest neighbors height differences, which have to be satisfied after the deposition. The 
randomly chosen deposition site is rejected (the atom is not deposited) if these restrictions are not satisfied. As a 
consequence, the number of deposition attempts does not coincide with the number of successful depositions in the 
KK model, although they are linearly related. 

All conserved growth models satisfy the conservation law 




(8) 




1-/3. 



(9) 



III. ATOMISTIC GROWTH MODELS 



dh(r,t) 
dt 



V-j(r,f)+77(r,t), 



(10) 



where j is the surface current and r\ is the noise term. Using different expressions, dictated primarily by symmetry 
considerations, for the current j, one can obtain all the conserved Langevin equations discussed in Sec. Ill Al The 
atomistic growth models considered in our work provide discrete realizations of these continuum growth equations. 
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It is known that some of the discrete growth models we study here have complicated transient behavior |42L |43| . 
For this reason, obtaining the dynamic scaling exponents that show the true universality classes of these models is 
often quite difficult. To make this task easier, the noise reduction technique 0,0] was introduced in simulations of 
such models. It has been shown |39| that this technique helps in suppressing high steps in the models and reduces 
the corrections in the scaling behavior, so that the true asymptotic universality classes of the growth models can be 
seen in simulations that cover a relatively short time. This makes it interesting to examine whether the persistence 
probabilities in these discrete models also exhibit similar transient behavior, and whether the noise reduction technique 
can help in bringing out the true persistence exponents of these models. To investigate this, we have applied the noise 
reduction technique to some of the discrete models studied in this paper. 

The noise reduction technique can be easily incorporated in the simulation of any discrete growth model by a small 
modification in the diffusion process [39l |. When an atom is dropped randomly, the regular diffusion rules for the 
growth model are applied and the final allocated site is chosen. Instead of adding the atom at that final site, a counter 
at that site is increased by one but the height of that site remains unchanged. When the counter of a lattice site 
increases to the value of a pre-determined noise reduction factor, denoted by m, the height at that lattice site is 
increased by one and the counter of that site is reset back to zero. The value of the noise reduction factor m should 
be chosen carefully. If m is too small, the suppression of the noise effect is not enough and the true universality class 
is not seen. However, is m is too large, the kinetically rough growth becomes layer-by-layer growth |46| and the 
universality class of the model cannot be determined. 

The atomistic models considered in our work are defined below. 

(i) Family model: The Family (F) model 2] is an extensively studied SOS discrete stochastic model, rigorously 
known to belong to the same dynamical universality class as the EW equation. It allows the adatom to explore within 
a fixed diffusion length to find the lattice site with the smallest height where it gets incorporated. If the diffusion 
length is one lattice constant (this is the value used in our simulations), the application of this deposition rule to a 
randomly selected site j involves finding the local minimum height value among the set: h(xj-x), h(xj) and h(xj + i) 
(in (l+l)-dimensions). The height of the site with the minimum height is then increased by one. 

(ii) Larger curvature model: The Kim-Das Sarma model p| is a more complex one which allows the atomic 
surface current j to be written as a gradient of a scalar field K, j = — VK, which can depend on h, V 2 /i, |V/i| 2 and 
so on. In the particular case when K = — V 2 /i, one obtains the so-called larger curvature (LC) model. As the name 
suggests, the diffusion rules applied to a randomly selected site j allow the adatom to get incorporated at the site in 
the neighborhood of site j where the local curvature (given by h(xj+i) + h(xj—i) — 2h(xj) in (l+l)-dimensions) has 
the largest value. The LC model asymptotically rigorously belongs to MH universality class described by Eq. 

(iii) Wolf— Villain model: The diffusion rules of the Wolf-Villain (WV) model 4] allow the adatom to diffuse 
to its neighboring sites in order to maximize its local coordination number which, for the (1+1 )-dimcnsional case, 
varies between 1 and 3 when the bond with the atom lying below the site under consideration is taken into account. 
In contrast to the F model, in this case the surface develops deep valleys with high steps almost perpendicular to 
the substrate. For the range of times and sample sizes used in the present study, the WV model may be considered 
to belong to the MBE universality class 0,03 described by Eq.JSJ. However, recent studies [3^. l47j have shown 
that the asymptotic universality class of this model in (l+l)-dimensions is the same as that of the EW equation. In 
contrast, in (2+l)-dimensions, studies based on the noise reduction technique have revealed that the WV model 
exhibits at very long times unstable (mounded) dynamic universality which cannot really be described by any of the 
continuum equations (Eqs. ©-(0) given above. 

(iv) Das Sarma Tamborenea model: The Das Sarma- Tamborenea (DT) model || is characterized by diffusion 
rules that are slightly different from those in the WV model. In this case, the diffusing atom tries to increase its 
coordination number, not necessarily to maximize it. For example, if a randomly selected deposition site has its local 
coordination number equal to 1 (i.e. no lateral neighbor in (l+l)-dimensions), and the two neighbors of this site have 
coordination numbers equal to 2 and 3, the deposited atom does not necessarily move to the neighboring site with 
the larger local coordination number: it moves to one of the two neighboring sites with equal probability (the atom 
would necessarily move to the site with coordination number 3 in the WV model). This minor change in the local 
diffusion rules actually changes the asymptotic universality class: the (l+l)-dimensional DT model belongs to the 
MBE universality class |39l l48| corresponding to the nonlinear continuum dynamical equation of Eq. However, 
the (2+l)-dimensional DT model asymptotically belongs to the EW universality a t very long times. 

(v) Controlled Kim— Das Sarma model: The Kim-Das Sarma model mentioned above provides a discrete 
realization of the continuum equation of Eq. (J5J if the scalar field K is chosen to be K = — V 2 h + \22iyh) 2 . However, 
the discrete treatment of the spatial gradients produces strong instabilities in the growth process due to uncontrolled 
growth of isolated structures, such as pillars or grooves. These instabilities can be easily controlled by introducing 
higher order nonlinear terms 0]. We call this new model the controlled Kim-Das Sarma (CKD) model. In this 
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model, the scalar field K is chosen to be K — — V h + A22/OV/1I ), where the nonlinear function / is given by 

1 _ p -c\Vh\ 2 

f(\Vh\ 2 ) = , (11) 

c 

with c > being the control parameter. The CKD diffusion rules for a randomly chosen deposition site, j, imply the 
minimization of the scalar field K , using the standard discretization scheme for the lattice derivatives V 2 /i and V/i: 

(V 2 % = h(x j+ i) + h(xj-i) - 2h(xj), (12) 



= -[h(x j+1 ) - hix^x)} 2 , (13) 

in (l+l)-dimensions. By carefully choosing the values for c and A22 [HI, one can remove the nonlinear growth 
instabilities completely and ensure an overall behavior of the CKD model similar to that of the DT model. 

(vi) Kim Kosterlitz and Kim— Park— Kim models: For completeness, we also present in this paper the results 
for the RSOS Kim-Kosterlitz (KK) and Kim-Park-Kim (KPK) § models which are known to belong asymp- 
totically to the KPZ and MBE universality classes, respectively. The common feature of these two models is the 
replacement of the usual diffusion rules of the SOS models described above by local restrictive conditions controlling 
nearest-neighbor height differences. 

In the KK model, deposition sites are randomly chosen, but the incorporation of the adatoms into the substrate is 
subject to a specific restriction: the deposition event occurs if and only if the absolute value of the height difference 
between the randomly selected deposition site, j, and each of its nearest-neighboring sites remains smaller than or 
equal to a positive integer n after deposition (our simulations were done for n = 1). If this strict constraint is not 
satisfied, the attempted deposition of an adatom is rejected, and the random selection of the deposition site is repeated 
until the deposition is successfully done. Since every attempt to deposit an adatom is not successful, the definition of 
"time" in this model is not quite the same as that in the other models where every deposition attempt leads to the 
incorporation of a new adatom in the growing film. In the KK model, the "time" is equivalent to the average height, 
which is not the same as the number of attempted depositions per site (these two quantities are the same in the other 
models considered here). The KK model is known to belong to the KPZ universality class, and in fact provides the 
most numerically efficient and accurate method for calculating the KPZ growth exponents. 

Kim et al. @ discovered that a slight change in the algorithm for choosing the incorporation site transforms the 
KK model into a new one, the KPK model, that belongs to the MBE universality class. The change consists of 
extending the search for appropriate incorporation sites (i.e, sites where the constraint on the absolute values of the 
nearest-neighbor height differences would be satisfied after the incorporation of an adatom) to the neighbors of the 
originally selected deposition site j. If the original site does not satisfy the constraint, then the neighboring sites 
(j ± 1 in (1+1) -dimensions) are checked, and an adatom is incorporated at one of these sites if the incorporation 
does not violate the constraint. Otherwise, the search is extended to the next-nearest-neighbors of j, and so on until 
a suitable incorporation site is found. We mention that in our implementation of this process, if, for example, both 
the sites j — k and j + k are found to be suitable for incorporation, then one of them is chosen randomly without 
any bias. Application of this algorithm in (2+l)-dimensions involves extending the search for suitable incorporation 
sites to those lying inside circles of increasing radii around the randomly selected deposition site j. The diffusion and 
incorporation rules of the KPK model @ lead essentially to a conserved version of the Kim-Kosterlitz RSOS model 
0, and as such the continuum growth equation corresponding to the KPK model is the conserved KPZ equation 
(with nonconserved noise), which is precisely the MBE equation; Eq. J5J is the conserved version of Eq. Q with 
nonconserved noise in both. 



IV. SIMULATION RESULTS AND DISCUSSION 



A. Persistence exponents in (1+1) dimensions 

Simulations for (l+l)-dimensional discrete growth models were carried out for (3 — 1/4, 3/8 and 1/3. The value 
j3 = 1/4 corresponds to the F model that has a relatively small equilibration time (of the order of L 2 ). The remaining 
conservative models, characterized by /3 = 3/8 (LC) and ~ 1/3 (WV, DT, CKD and KPK), have a much slower 
dynamics (with z values 4 or 3). So their corresponding equilibration time intervals, required for the interface 
roughness to reach saturation, are of the order of L 4 and L 3 , respectively. For this reason, the largest values of L for 
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Growth model 


L 


el 


el 


P 


Universality class 


F 


10 6 


1.57 ±0.10 


1.49 ±0.10 


0.25 ±0.01 


EW 


KK 


5 x 10 4 


1.68 ± 0.02 


1.21 ±0.02 


0.33 ±0.01 


KPZ 


LC 


10 4 


0.84 ± 0.02 


0.84 ± 0.02 


0.37 ±0.01 


MH 


WV 


10 4 


0.94 ± 0.02 


0.98 ±0.02 


0.37 ±0.01 


MBE 


DT 


10 4 


0.95 ± 0.02 


0.98 ±0.02 


0.38 ±0.01 


MBE 


CKD 


10 4 


0.98 ± 0.02 


0.93 ±0.02 


0.35 ±0.01 


MBE 


KPK 


10 4 


1.04 ±0.02 


1.01 ±0.02 


0.31 ±0.01 


MBE 



TABLE I: Positive and negative persistence exponents, 8+ and 0-, for the transient (T) regime, measured for seven different 
discrete growth models (identified in the first column) using kinetic Monte Carlo simulations with relatively large system 
sizes (L). The measured growth exponent, /3, and the universality class of the model are indicated in the last two columns, 
respectively. 



Growth model 


L 


el 


8 T _ 


9% 


8 S _ 


P 


F 


10 3 


1.67 ±0.10 


1.56 ±0.10 


0.78 ±0.02 


0.76 ± 0.02 


0.25 ±0.01 


KK 


5 x 10 2 


1.70 ±0.02 


1.27 ±0.02 


0.71 ±0.02 


0.71 ± 0.02 


0.30 ±0.01 


LC 


40 


0.98 ±0.02 


0.96 ±0.02 


0.67 ±0.02 


0.67 ±0.02 


0.32 ±0.01 


WV 


40 


0.94 ±0.02 


0.99 ±0.02 


0.65 ± 0.02 


0.70 ±0.02 


0.35 ±0.01 


DT 


40 


0.98 ±0.02 


1.01 ±0.02 


0.64 ±0.02 


0.72 ±0.02 


0.36 ±0.01 


CKD 


40 


1.11 ±0.02 


0.99 ±0.02 


0.78 ± 0.02 


0.66 ±0.02 


0.33 ±0.01 


KPK 


2 x 10 2 


1.16 ±0.02 


1.09 ±0.02 


0.70 ±0.02 


0.68 ±0.02 


0.28 ±0.01 



TABLE II: Positive and negative persistence exponents, 8+ and 8-, for the transient (T) and the steady state (S) regimes of 
our seven different discrete growth models, obtained from simulations with relatively small samples sizes (L). To illustrate the 
effects of reduced system sizes on the measured exponents, we have shown the values of /3 obtained from these simulations in 
the last column. 

which the steady state could be reached in reasonable simulation time are considerably shorter in these models than 
in the F model. The fastest equilibration occurs in the KK model (j3 = 1/3) where z = 3/2. 

In calculations of the transient persistence probabilities, the initial configuration of the height variables is taken to 
be perfectly flat, i.e. hj(to) = (j = 1, L). The lattice size was in the range 10 4 < L < 10 6 , and the duration of the 
deposition process, measured in units of number of grown monolayers (ML), was ~ 10 3 . The results were averaged 
over ~ 10 3 independent runs. For measurements in the steady-state situation, a saturation of the interface roughness 
was first obtained by depositing a large number (of the order of L z ) of monolayers and subsequent time evolution 
from one of the steady-state configurations obtained this way was used for measuring the persistence probabilities. 
A much smaller lattice length (L = 1000 for the F model, L = 500 for the KK model, L = 200 for the KPK model, 
and L — 40 for the LC, WV, DT and CKD models) was used in these calculations in order to reach the steady-state 
saturation within reasonable simulation times. 

The positive (negative) persistence probabilities in both transient and steady-state regimes were obtained as the 
fraction of sites that maintain the values of their heights persistently above (below) their initial values, averaged over 
a large number (~ 10 4 ) of independent runs. The persistence exponents were obtained from power law fits to the 
decay of these probabilities, as shown in Figs. ^0] and HJHE1 for the transient and steady-state regimes, respectively. 

For all the models studied here, we have also measured the value of the growth exponent f3 in both transient 
and steady-state simulations. Since the latter simulations were carried out for smaller values of the system size 
L, these measurements provide useful information about the dependence of the measured exponent values on the 
lattice size. Similar information is also provided by the values of the transient persistence exponents obtained from 
measurements in the initial stage of the steady-state simulations. The transient exponent values obtained from 
the large-L simulations are listed in Table [IJ and both transient and steady-state exponent values obtained from 
simulations of relatively small samples are shown in Table ITT1 The measured values of the growth exponent j3 are also 
shown in these Tables. 

Estimation of the probable error in the measured values of the growth and persistence exponents is a delicate task 
(and surely depends on precisely how the exponent error is defined) , since there is not a traditional accepted method 
to evaluate the error in dynamical simulations. To solve this problem we did the following simulations. We decreased 
the number of independent runs used for the averaging procedure by a factor of 2, keeping the size of the system 
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FIG. 1: Transient persistence probability for the (l+l)-dimensional linear F and LC growth models. As expected, the positive 
and negative persistence probabilities are identical in these models. The system size is L = 10 6 for the F model and L = 10 4 
for the LC model, and an average over 10 3 independent runs was performed. The slopes of the double-log plots yield the values 
of the transient persistence exponents shown in Table Q 



constant. Under these circumstances, we have measured the exponents corresponding to the two different numbers of 
independent runs and the differences between the obtained values of the exponents were used as error estimates for j3 
and 9, respectively. Approximately the same size of the error bar was obtained from the estimations of fluctuations 
in the value of the local slope of the double-log plots. We have also noticed that a reduction of the lattice size 
(imposed for the steady-state persistence calculations) produces lower values of the growth exponents, as shown in 
Table ITT1 This is because the downward bending (approach to saturation) of double-log width versus time plots occurs 
at shorter times in simulations of smaller systems. However, the smaller-L simulations seem to lower the measured 
values of the growth exponents by a maximum of about 10% percent. So we conclude that this effect is not dramatic 
and that the steady-state results reported below are reliable. 

The measured values of (3 agree reasonably well with the expected ones (see Section III A|) within their errors. As 
expected, the agreement is better in the case of larger values of L. For the larger-!/ simulations (L ~ 10 4 ), we have 
found that the growth exponents of the F, LC and KK models are in excellent agreement with their corresponding 
expected values of 1/4, 3/8 and 1/3, respectively (see Tabled . The DT and WV models are found to behave similarly 
at early (transient) stages of their interface growth, at least in (l+l)-dimensions, their growth exponents being: 
fiwv ~ 0.37 and (3dt ~ 0.38. The closeness of these values to the value of 3/8, which corresponds to the MH 
universality class, suggests that the nonlinear term that appears in the associated dynamic equation (i.e. Eq. JSJ) has 
a very weak effect for the range of lattice sizes used in our study. In addition, we have found that the CKD model 
characterized by the nonlinear coefficient A22 = 2 and control parameter c = 0.02 has a growth exponent (3ckd ~ 0.35, 
in agreement with Ref. |14|. These particular parameter values ensured the elimination of any interfacial instability, 
thus allowing a calculation of the steady-state persistence properties. Regarding the conserved KPK model, we have 
observed that the growth exponent has a value that is slightly smaller than 1/3, a result that agrees with Ref. 

The temporal behavior of the transient persistence probability in our models is shown in Figs. EE! From these 
measurements, we obtained the transient persistence exponents by fitting the linear middle regions (excluding the 
small-f and large-i ends, typically using the data for 20 < t < 800) of the double-log plots to straight lines. As 
expected, due to the invariance of the interfaces of the F and LC models (which are characterized by linear continuum 
equations) under a change of sign of the height variables, we obtained equal positive and negative transient persistence 
exponents within the error bars, as displayed in Fig. ^ However, we mention that the F model has a rather slow 
convergence of the positive and negative exponents towards their long-time value of ~ 1.55 observed in much longer 
simulations. The results for F and LC models, that correspond to (3 = 1/4 and 3/8, respectively, agree well with the 
values reported by Krug et al. [23|. The same level of agreement is also found in the case of the KK model [H, 
shown in Fig. [3 for which the transient persistence exponents are 9^ « 1.68 and « 1.21 in (l+l)-dimensions. We 
note that the negative persistence probability has a slower decay than the positive one. This is due to the constant 
coefficient, A2, of the nonlinear term |V/i(r, t)\ 2 of the KPZ equation (which provides a continuum description of the 
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FIG. 2: Positive and negative transient (bottom two curves) and steady-state (top two curves, mostly overlapped) persistence 
probabilities for the (l+l)-dimensional RSOS KK model. The faster decay of the positive persistence probability in the 
transient regime is due to the negative sign of A2 in the equivalent continuum equation of Eq. JHJl. In the transient case, 
systems of size L = 5 x 10 4 were averaged over 5 x 10 3 independent runs. The steady-state simulation was done for L = 500 
and a similar average was performed. 

KK model) having a negative sign p^j . 

For the models described by the fourth-order nonlinear MBE equation (i.e. WV, DT, CKD and KPK models), we 
expect to find different positive and negative transient persistence exponents due to the fact that their morphologies 
violate the up-down interfacial symmetry with respect to the average level. No information about how different 
these two exponents should be is available in the literature. In most of these growth models, we observe that the 
two exponents are not very different from each other, especially during the transient regime. Fig. [3] shows the 
transient regime results for DT and WV models, which are indeed very similar - their persistence probability curves 
have almost identical behavior. We note here that the negative persistence probability has a faster decay than the 
positive persistence probability. This indicates a negative sign of A22, the coefficient that multiplies the nonlinear 
term V 2 \Vh(r, t)\ 2 of the MBE equation. However, the relative order of the values of these exponents is reversed 
when A22 > 0, which is the case in the CKD and KPK models, as shown in Fig.0] To clarify this aspect, we show in 
Fig- El the interfacial morphologies of DT and CKD models. We used a lattice of L — 10 4 sites (but only a portion of 
1000 sites is shown in each case) and the displayed configurations correspond to a time of 10 3 ML. The interface of 
the DT model is characterized by deep grooves, while the profile in the CKD model exhibits the distinct feature of 
high pillars. Both morphologies display strong up-down interfacial asymmetry, but their representative features (i.e. 
deep grooves and high pillars) are opposite in "sign" , indicating a reversal of the sign of the coefficient A22 (note that 
a reversal of the sign of A22 in Eq. (JjjJ is equivalent to changing the sign of the height variable h(r, t)). 

As summarized in Table U the DT, WV and CKD models show very similar values for the transient persistence 
exponents when the above mentioned effect of the sign of A22 is taken into account. However, some deviation from 
the exponent values for this group of models is observed in the RSOS KPK model which shows the smallest difference 
between the positive and negative persistence exponents. Finite size effects appear to be stronger in this case. These 
effects also cause an increase in the measured values of the persistence exponents above the expected values. A similar 
behavior is found in the steady-state results as well, as described below. 

Our calculations of WV, DT, CKD and KPK persistence exponents illustrate the feasibility of studying this type 
of nonequilibrium statistical probabilities for a large class of nonequilibrium applications described by nonlinear 
dynamical equations. Until now, the only nonlinear equation for which persistence exponents have been calculated 
|24j is the KPZ equation which is arguably the simplest nonlinear Langevin equation. Further, the nonlinearity in 
the KPZ equation becomes irrelevant in the steady-state regime in (14T) -dimensions. So, the effects of nonlinearity 
are not reflected in the steady-state persistence behavior of (l+l)-dimensional KPZ systems. An immediate concern 
would be that more complex nonlinear dynamic equations might be less approachable from the point of view of 
persistence probability calculation. Our results for four nonlinear models eliminate this possibility and illustrate the 
applicability and usefulness of persistence probability calculations in the study of surface fluctuations. 
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FIG. 3: Positive and negative transient persistence probabilities for the (1+1) -dimensional nonlinear DT and WV growth 
models. We note that despite the difference in their local diffusion rules, these two models behave identically as far as the 
transient persistence probability is concerned. The curves corresponding to the DT model have been shifted upward in order 
to avoid a complete overlap of the plots for the two models. The system size is L = 10 4 and an average over 10 3 independent 
runs was performed. The slopes of the double-log plots yield the transient persistence exponents given in Table Q] 



Figures EHH1 display our results for the steady-state persistence probabilities. The values of the growth and persis- 
tence exponents obtained from the steady-state simulations are summarized in TableHTl The values of the steady-state 
persistence exponents in the F and LC models, corresponding, respectively, to the V 2 and V linear equations, (see 
Fig. EJ) are consistent with the values of the corresponding growth exponents (as predicted by Eq. @) obtained from 
the same small-!/ simulations. For the WV and DT models, as shown in Fig.[7| we obtain very similar positive and 
negative persistence exponents. In the case of the KK model wc find, as expected, identical positive and negative 
exponents (9± « 0.71), as shown in Fig. [21 

Among the models belonging to the MBE universality class, the KPK model exhibits steady-state persistence 
exponents that are systematically higher than the ones obtained for the remaining three (WV, DT and CKD) models. 
Our study of the dynamical scaling behavior of the KPK model indicates that both a (~ 0.9) and z (~ 2.9) in this 
model are reasonably close to the expected values, in agreement with Ref. |8|. Therefore, the reason for the differences 
between the values of the steady-state persistence exponents for the KPK model and those in the other models in 
the MBE universality class is unclear. This discrepancy may very well be arising from subtle differences in finite 
size (and time) effects in the simulations for persistence exponents and dynamic scaling. Further investigation of the 
applicability of this RSOS model in understanding MBE growth is beyond the scope of the present study. 

As shown in Fig. |S] a nice illustration of the presence of nonlinearity in the underlying dynamical equation is 
provided by the steady-state persistence exponents of the CKD model, characterized by distinct values, w 0.78 
and Q s _ » 0.66, of the positive and negative exponents. Although one must take into account the fact that these 
values might be slightly overestimated (by approximatively 5%) due to the smallness of the sample sizes used in the 
steady-state simulations, these exponents provide a good qualitative account of the nontrivial up-down asymmetric 
persistence behavior expected for nonlinear models belonging to the MBE universality class. 

Next we investigate the influence of small sample sizes on the measured values of the persistence exponents. We 
mainly used the DT model to answer this question and we pursued the following two tests. First we decreased the 
size of the system from 10 4 to 100 and then to 40 and found the values of the growth and persistence exponents, as 
summarized in Table ITTT1 We note that as the lattice length decreases to 40, the persistence exponents increase by ~ 
2%, while the growth exponents increase by ~ 5%. As a second test, we have applied the noise reduction technique to 
both the DT and WV models. It has been shown |39j that a noise reduction factor of m = 5 helps the DT model to 
recover quite accurately the universal exponents corresponding to the MBE universality class. In addition, the noise 
reduced WV model exhibits, at late evolution times, its true EW asymptotic universality, which is difficult to observe 




FIG. 4: Positive and negative transient persistence probabilities for the (l+l)-dimensional CKD (the upper two curves) and 
RSOS (the lower two, almost overlapped curves) models that belong to MBE universality class. In both cases the system size 
was L = 10 4 and an average over 10 3 independent runs was performed. The slopes of the double-log plots yield the transient 
persistence exponents given in Table 
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0.98 ± 0.02 
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TABLE III: Transient positive and negative persistence exponents, obtained for the DT model with different system sizes 
(I/). The effect of the system size on the measured growth exponent, /3, is displayed in the last column. No result for steady- 
state persistence exponents is available for system sizes larger than ~ 100, due to the impossibility of reaching saturation of 
the interface width for such values of L in time scales accessible in simulations. The results shown here were averaged over 500 
(for L = 10 4 ), 5 x 10 4 (for L = 100) and 10 5 (for L = 40) independent runs. 



without applying noise reduction. Therefore, the DT model with the appropriate noise reduction factor is expected 
to provide the correct persistence exponents associated with the fourth- order nonlinear dynamical equation for MBE 
growth. The results obtained from the simulations with noise reduction are summarized in Table ITVl We notice that 
the noise reduction scheme produces only a minor change in the persistence exponents and in addition, the results 
obtained for m = 5 agree within the error bars with those for the CKD model. We, therefore, conclude that the 
noise reduced DT model and the discrete CKD model provide a good representation of the MBE universality class, 
characterized by two different steady-state persistence exponents: 6+ ~ 0.66 (positive persistence) and ~ 0.78 
(negative persistence). These nontrivial persistence exponents for this class have not been reported earlier, and it 
would be useful to check these results from further theoretical or experimental studies. Regarding the noise reduced 
WV model we mention that the convergence of 9 s towards the expected value of 3/4 is rather slow in the case of the 
positive exponent and probably a higher value of the noise reduction factor would be necessary to reveal the true EW 
universality. We did not explore this technical issue any further. 

We note that among the positive and negative steady-state persistence exponents for these nonlinear growth models, 
the smaller one (for example, the positive exponent in the DT model or the negative exponent of the CKD model), 
turns out to be close to (1 — 0). In the next subsection, we show analytically that this relation between the smaller 
steady-state persistence exponent and the dynamic growth exponent is, in fact, exact. Our numerical studies suggest 
a connection of this result with the morphology that develops in the steady-state regime. As shown in Fig. [SJ the 
characteristic feature of the DT morphology is the presence of deep grooves, while the CKD model exhibits high 
pillars. Loosely speaking, in the case of the DT model, we expect the relation of Eq. @ to be more likely to be 
satisfied by the positive persistence exponent than the negative one because the preponderant grooves, responsible for 
the negative persistence exponent, represent the effects of the nonlinearity of the underlying MBE dynamics. More 
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FIG. 5: Morphologies of the (l+l)-dimensional DT (top) and CKD (bottom) stochastic models for L = 10 4 (only a portion 
of 1000 sites is shown) and t = 10 3 ML. In the DT model, we notice a breaking of up-down symmetry due to the formation of 
deep grooves, while in the CKD model, the representative asymmetric feature corresponds to high pillars. 

work is clearly needed for a better understanding of the possible relationship between such "nonlinear" features of 
the interface morphology and the value of the persistence exponent. 

B. An exact relation between steady-state persistence exponents and the growth exponent 



As mentioned earlier, for interface heights h(r,t) evolving via a Langevin equation that preserves (h — * —h) 
symmetry (for example, any linear Langevin equation), the steady state persistence exponents satisfy the scaling 
relation 0+ = — 1 — 0, where j3 is the growth exponent |2^| . In this subsection, we derive a generalized scaling 
relation, 



which is valid even in the absence of (h 
known result pl|. 6>f = 6 s = 1 - 0. 



(3 = max [l - 6^,1 - s ] , (14) 
-h) symmetry. When this symmetry is restored, Eq. (|14*|) reduces to the 
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FIG. 6: Positive and negative steady-state persistence probabilities for the (l+l)-dimensional F and LC models which are 
governed by linear continuum dynamical equation. The temporal decay of the persistence probability is slower in the LC model 
which has a larger growth exponent ({3lc = 3/8, /3p = 1/4). We used L = 1000 and t — 4 x 10 6 ML for the F model, and 
L — 40, to — 10 6 ML for the LC model. The displayed results were averaged over 5000 independent runs. The measured slopes 
of the double-log plots yield the steady-state persistence exponents shown in Table ITTI 
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FIG. 7: Steady-state persistence probabilities for two (l+l)-dimensional models in the MBE universality class - the DT and 
WV models. As in the transient case, these two models exhibit almost identical persistence behavior in the steady state. The 
effects of the nonlinearity in their continuum dynamical description are not very prominent for the small lattice sizes considered 
here. For the data shown, systems of size L = 40 were equilibrated for to = 10 5 ML, and the results were averaged over 5000 
independent runs. The persistence plots for the DT model have been shifted up in order to make them distinguishable from 
the WV plots. The measured slopes of the double-log plots yield the steady-state persistence exponents shown in Table ITTI 




FIG. 8: Double-log plots of the steady-state persistence probabilities of (1+1) -dimensional MBE class CKD and KPK (shifted 
up by a constant amount) models. While the KPK model does not show a clear effect of nonlinearity in the values of the 
persistence exponents, the CKD model shows positive and negative persistence exponents that are clearly different from each 
other. Systems of size L — 40 (CKD) and L — 200 (KPK) were equilibrated for to ~ 10 5 ML. The results were averaged 
over 10 4 independent runs. The measured slopes of the double-log plots yield the steady-state persistence exponents shown in 
Table |II| 



Growth model 


m 


Of 


e s _ 


DT 


1 


0.64 ±0.02 


0.72 ±0.01 


DT 


5 


0.65 ±0.02 


0.77 ±0.01 


WV 


1 


0.65 ±0.02 


0.70 ±0.01 


wv 


5 


0.68 ±0.02 


0.75 ±0.01 



TABLE IV: Positive and negative persistence exponents, Q±, for the steady state of the DT and WV models for two different 
values of the noise reduction factor, m. Systems of size L = 40 were equilibrated for 10 5 ML and the results were averaged 
over 5000 independent runs. 



To derive the relation in Eq. (|14fl . we start with a generic interface described by a height field h(r, t) and define the 
relative height, u(r,t) = h(r,t) — h(r,t) where h(r,t) = J h(r,t)dr/V is the spatially averaged height and V is the 
volume of the sample. Let us also define the incremental auto-correlation function in the stationary state, 

C{t,t')= lim ([u(r,t + t )-u(r,t' + t )] 2 }. (15) 

to — >oo 

It turns out that for generic self-affine interfaces (which do not have to be Gaussian), this function C(t, t') depends only 
on the time difference \t — 1'\ (and not on the individual times t and t') in a power-law fashion for large \t — 1'\ p3.l49|. 



C{t,t')~\t-t'\^, (16) 

where (3 is the growth exponent. 

This particular behavior of the auto-correlation function in Eq. Hlrj|) is typical of a fractional Brownian motion 
(fBm). A stochastic process x(t) with zero mean is called an fBm if its incremental correlation function C(t-\ , fe) = 
{[x(t\) — xfo)} 2 ) depends only on the time difference \t\ — *2 1 in a power-law fashion for large arguments 

C(ti, h) = ([x(ti) - x{t 2 )] 2 ) ~ \t x - t 2 \ 2H , (17) 



where < H < 1 is called the Hurst exponent of the fBm. For example, an ordinary Brownian motion which evolves 
as dx/dt = rj{t) where r](t) is a Gaussian white noise with zero mean and a delta function correlator, satisfies Eq. i|17|) 
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with H — 1/2. Thus an ordinary Brownian motion is a fBm with H = 1/2. It follows clearly by comparing Eqs. i|16|) 
and H17|) that the relative height u(r, t) of a generic interface at a fixed point r in space, in its stationary state, is also 
a fBm with Hurst exponent, H = (3. Note that an fBm is not necessarily Gaussian. 

We are then interested in the 'no return probability' to the initial value of the fBm process u(r, t). So, the relevant 
random process is Y(r, t) = u(r, t + to)— u(r, to). Clearly, Y(r, t) is also a fBm with the same Hurst exponent /3 since 
the incremental correlation function of Y is the same as that of u(r, t). We are then interested in the zero crossing 
properties of the fBm Y(r, t). Now, consider the process Y(r,t) as a function of time, at a fixed point r in space, 
from time to to time to + 1 where to — > oo. There are two types of intervals between successive zero crossings in time, 
the '+' type (where the process lies above 0) and the ' — ' type (where the process lies below 0). 

In general, the statistics of the two types of intervals are different. Only, in special cases, where one has the 
additional knowledge that the process Y(r,i) is symmetric around (i.e., processes which preserve the (h — » —h) 
symmetry), the + and — intervals will have the same statistics. For such cases, a simple scaling argument was given 

in Ref. [23| to show that the length of an interval of either type has a power-law distribution, Q(t) ~ r -1-9 (for 
large r) with 9 s = 1 — H =1 — (3. Note that this relation between the persistence exponent and the Hurst exponent is 
very general and holds for any symmetric fBm, i.e., any stochastic process with zero mean (not necessarily Gaussian) 
satisfying Eq. (|17fl . Recently, other applications of this result have been found [EH IH^ . For general nonsymmetric 

processes, however, one would expect that Q+(t) ~ t~ 1 ~ 9± for large r, where 0+ and 0*L are, in general, different. 
Here we generalize this scaling argument of Ref. |23j (derived for a symmetric process) to include the nonsymmetric 
cases and derive the result in Eq. I|14l) . 

The derivation of Eq. (jl4|l follows more or less the same line of arguments as that used in Ref. |2^| for the symmetric 
case. Let P(Y, r) denote the probability that the process has value Y at time r, given that it starts from its initial 
value at r = 0. Then, it is natural to assume that the normalized probability distribution P(Y, r) has a scaling 
form, 

P(Y,T) = ^-f(^-), (18) 

where cr(r) is the typical width of the process, a 2 (t) = (Y 2 {t)). It follows from Eq. JE) that er(r) ~ t 13 for large r. 
The scaling function f(z) is a constant at z = 0, /(0) ~ 0(1) (note that, in general, f(z) is not a symmetric function 
of z) and should decrease to as z — > ±oo. So, given that a zero occurs initially, the probability p(r) — P(0,r) that 
the process will return to after time r (not necessarily for the first time) scales as 

P {r)~±- r r-e, (19) 

as r — > oo. This function p(r) indeed is the density of zero crossings between r and r + dr. Thus, the total number 
of zeros upto a time T is simply the integral, 

N(T)= f p(T)dr~T 1 -' 3 , (20) 
Jo 

for large T. 

Next, we relate the persistence probabilities to the number of zeros. Let P±(t) denote the probabilities that the 
process stays positive (or negative) over the interval [0, r], given that it started from a zero. By definition, we have 
P±{t) ~ t~ 6± for large r. Then, Q±(r) = —dP±(r)/dT ~ r _1_e ± (as r — > oo) denote the probabilities that the 
process will cross zero next time (from the positive or the negative side respectively) between time r and r + dr. 
Thus, Q±(t) are also the distribution of intervals of the two types of length r. 

Now, consider a total length of time T. Let N(T) denote the total number of intervals in this period, half of them 
are + types and the other half — types, N±(T) = N(T)/2. Let, n±(r) denote the number of ± intervals of length 
t within the period T. Thus, the fraction of + (or — ) intervals of length r, ti±(t)/N±(t), by definition, are the two 
distributions Q±(t) provided T is large. Thus, for large T, we have 

n±(r,T) = ^pQ ± (r) ~ ^V(T)r~ 1 ~ 9 ±, (21) 

for 1 << t < T. On the other hand, we have the length conservation condition (the total length covered by the 
intervals must be T), 



/ dr t [n + (t) + ?i_ (t)] = T. 
Jo 



(22) 
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Substituting the asymptotic behavior of n± (r) in Eq. (|21|l into the left-hand side of Eq. (1221 , we get 



N(T) 



T i- 



T 



1 - i 



1 - 



oc T. 



(23) 



We next use N(T) ~ T 1_/3 for large T from Ea.pty. This gives, for large T, 



(24) 



Taking T — > oo limit and matching the leading power of T in Eq. (|24|l . we arrive at our main result in Eq. I] 14(1. Note 
that in the above derivation we have implicitly assumed a small-i cut-off and focused only on the distribution of large 
intervals. Our numerical results obtained for a class of nonlinear interfaces in both (1 + 1) and (2 + 1) dimensions 
(see Sec. IIV Dl below') are consistent with the analytical result in Eq. (|14fl . 



C. Dependence of persistence probabilities on the initial configuration 

We present in this section some surprising simulation results about the dependence of the persistence behavior 
(specifically, the values of the persistence exponents) on the choice of the initial configuration. In particular, we 
show that the steady-state exponents may be obtained with a fair degree of accuracy from simulations in which the 
interface has not yet reached the steady state. We also present some results that have bearing on the measurability 
of the transient persistence exponents from experimental data. 

We recall that in Section IIV Al the transient persistence exponents were measured from simulations in which the 
initial configuration was completely flat, corresponding to to — 0. To examine the dependence of the persistence 
probabilities on the choice of to, we evolved samples governed by F and DT atomistic diffusion rules for to = 10, 100 
and 1000 ML, starting from perfectly flat initial states and used the resulting configurations as starting points for 
measuring the persistence probability (the probability of the height at a given site not returning to its initial value 
at time to) as a function of t. We show the results of these simulations in Figs . 151 and 1 1 01 for the F and DT models, 
respectively. We find that even for the small value of i = 10 ML [see panel a)], the observed persistence probabilities 
do not exhibit power law decay in time with the transient persistence exponents, despite the fact that the expected 
condition |23| for transient behavior, t 3> t , is well satisfied in a large part of the range of t used in these simulations. 
These results point out a practical difficulty in obtaining experimental evidence for transient persistence behavior. 
Since perfectly flat initial configurations can hardly be achieved experimentally and experimental measurements are 
always started from a relatively rough substrate, the transient persistence exponents may very well not be measurable 
from experiments if the only way of measuring these exponents is to start from a perfectly flat morphology. 

As the value of to is increased to 100 ML, the persistence probabilities tend to show the expected power law behavior, 
as shown in Fig. EJb) for the F model and in Fig. HOf b') for the DT model. Most surprisingly, as shown in Figs.[!|fc) 
and llOf c'l. we find that for to — 1000 ML, one recovers precisely the power law behavior, P(to, to + t) oc t~ 6 , and the 
exponents are essentially the same as the previously obtained steady-state ones shown in Table[n] This investigation, 
thus, reveals the fact that a measurement of the steady-state persistence exponents does not require the preparation 
of an initial state in the long-time steady-state regime where the interface width has reached saturation: an initial 
state in the pre-asymptotic growth regime where the interface width is still increasing as a power law in time [as 
illustrated in Figs.|3d) and[TUfd)] is sufficient for measurements of the steady-state persistence exponents. A similar 
result was reported in Ref. |23| , but it was argued there that the measurement time t must be much smaller than to for 
steady-state persistence behavior to be observed. Our results show that the steady-state persistence exponents are 
found even if t is of the order of (or even slightly larger than) the initial time to- This observation has an important 
practical benefit: it implies that one can easily obtain accurate estimates of the steady-state persistence exponents 
using rather large systems (L ~ 10 4 ), and growing approximately up to t ~ 10 3 ML, instead of having to use the 
very large values (of the order of L z ) of to necessary for obtaining saturation of the interface width. At the same time, 
this observation also illustrates the above mentioned difficulty in obtaining the transient exponents from experimental 
measurements. 

To investigate the effects of random imperfections in the initial substrate (which are always present in experimental 
studies) on the persistence behavior, we carried out simulations in which particles were deposited randomly on a 
perfectly flat substrate for 10 ML and the resulting configuration was used for further depositions using the diffusion 
rules of the F and DT models. Persistence probabilities were calculated starting from the configurations obtained 
after the random deposition of 10 ML. Figures ITU and IT21 show the results for the F and DT models, respectively. We 
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FIG. 9: Log-log plots of the positive and negative persistence probabilities [panels (a)-(c)] for the F model, obtained using 
different values of the initial time to- Systems with L — 10 4 sites have been averaged over 500 independent runs. Persistence 
probabilities are computed starting from the configuration corresponding to: a) to = 10 ML. We do not find a clear power law 
decay of the persistence curves, b) to = 100 ML. As to increases, a clearer power law behavior is observed, c) to = 1000 ML. 
The power law decays are recovered and characterized by exponents in agreement with those corresponding to the steady-state 
regime: 0± w 0.75. d) Log-log plot of the interface width W as a function of time t (in units of ML). The value of the slope 
(equal to the growth exponent 0) agrees with the expected value, j3 = 0.25. 



find that even when the persistence calculation starts from a configuration characterized by random deposition, there 
is an indication that one can still obtain the steady-state exponents during the last decade of t where the growth 
exponent reaches the values characteristic of the diffusion rules of the specific (F or DT) model being considered. 
Indeed, in the time region where the growth exponents are (3 = 0.25 for the F model [see Fig. Illf b)] and (3 = 0.375 
for the DT model [see Fig. I12f b)] , we have calculated the persistence exponents and recovered values very close to the 
steady-state ones. These observations confirm our earlier conclusions about the relatively easy measurability of the 
steady-state persistence exponents and the difficulty in measuring the transient exponents in experimental situations. 



D. Persistence exponents in (2+1) dimensions 

Our calculations in (2+l)-dimensions make use of our observation (discussed above) concerning the possibility of 
obtaining the correct steady-state exponents from simulations that avoid the time consuming process of reaching 
the true steady state where the interface width has saturated. The result that the persistence exponents obtained 
from (l+l)-dimensional simulations using fairly small values of to and i ~ to are quite close to the steady-state 
values allows us to extract numerically the steady-state persistence exponents in (2+l)-dimensions using systems 
with reasonably large sizes. If one had to run systems of size L ~ 100 x 100 all the way to saturation in order to 
measure the steady-state persistence exponents, it would have been impossible to do the calculations within reasonable 
simulation time. In addition, decreasing the system size is not an acceptable solution because the results then become 
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FIG. 10: Positive and negative persistence probabilities [panels (a)-(c)] for the DT model, obtained using different values of 
the initial time to. Persistence probabilities are computed starting from the configuration corresponding to: a) to — 10 ML. 
As in the case of the F model, we do not find a clear power law for the persistence curves, b) to = 100 ML; c) to = 1000 ML: 
Power law decays are recovered and characterized by exponents that are approximately equal to those corresponding to the 
steady-state regime, 8+ ~ 0.64 and 9'1 ~ 0.71. d) Log-log plot of the interface width W as a function of time t (in units of 
ML). The slope gives a growth exponent of (3dt — 0.375. 

dominated by finite size effects. 

Simulations for (2+l)-dimensional discrete growth models were carried out for the F model (J3 = 0) and the DT 
model (/3 ~ 1/5). Simulations using systems of size L = 200 x 200 revealed that the growth exponents, obtained from 
averages over 200 independent runs, are (3 = 0.04 ± 0.01 and 0.20 ± 0.01 for the F and DT models, respectively, in 
agreement with Ref. [43. In the DT model, we noticed a crossover from the initial value of 0.26 to the asymptotic 
expected value of 0.20, indicating that no additional noise reduction technique is necessary for obtaining results that 
reflect the correct universality class of this model. For both F and DT models we calculate the transient and steady- 
state persistence probabilities by recording the fraction of sites which do not return to their initial height up to time 
t, as in the (l+l)-dimensional case. We used to = (perfectly flat initial state) in the calculation of the transient 
persistence probabilities, and three different values, such as to = 20 ML, 200 ML and 2000 ML for the F model, in 
the calculation of the steady-state exponents. 

We report the results for the transient probabilities just for the sake of completeness: the rapid decay of the 
persistence probabilities prevents us from obtaining accurate values of the associated persistence exponents. This fast 
decay of the transient persistence probability is a consequence of the reduced roughness of these higher dimensional 
models. This effect is particularly pronounced for the F model for which the persistence exponent is found to be 
larger than 6 and the persistence probability decreases rapidly to zero for any deposition time larger than ~ 60 ML, 
as shown in Fig. ^] We also observe that the transient values of the positive and negative persistence exponents 
in the DT model are roughly 3 times larger than the values obtained in the (1 + l)-dimensional case. The relative 
difference between the positive and negative persistence exponents remains approximately the same as that in the 
(l+l)-dimensional model. Our results for these (2+l)-dimensional persistence exponents are summarized in Table 
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FIG. 11: a) Positive and negative persistence probabilities for the F model. During the deposition of the first 10 ML, the 
growth process is random deposition. The diffusion rules of the F model are then used to evolve the interface. Persistence 
probabilities are computed starting from the configuration obtained after the random deposition of 10 ML. The positive and 
negative persistence exponents in the last growth decade are in the range 0.6 to 0.7, depending on the fitting region, b) 
Log-log plot of the interface width W as a function of t (in ML). The slope in the first decade of t is precisely the random- 
deposition value, /3 = 0.5. The second decade shows a crossover region where the systems undergoes a transformation towards 
a morphology governed by the F model diffusion rules, and the last decade is characterized by the expected growth exponent 
of the F model, p = 0.25. 
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FIG. 12: Positive and negative persistence probability curves for the DT model. During the deposition of the first 10 ML, 
the growth process is random deposition. The diffusion rules of the DT model are used to evolve the interface subsequently. 
Persistence probabilities are computed starting from the configuration obtained after the deposition of the first 10 ML. The 
positive and negative persistence exponents in the last growth decade are approximatively equal to 0.66 and 0.79 respectively, 
b) Log-log plot of the interface width W as a function of time t (in ML). Beyond the crossover regime, the last decade in t is 
characterized by the expected growth exponent of the DT model, p ~ 0.375. 
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We now focus on the steady-state persistence exponents which, as discussed above, are found using relatively small 
values of to and t as to- In Fig. ll4f a). we show that for to -C t (e.g. for to = 20 ML), the persistence probability of the 
F model does not exhibit a clear power law decay. However panels b) and c) of Fig. ^] reveal that once to becomes of 
the order of the measurement time t, the expected power law is recovered and in addition the steady-state exponent 
for the linear F model, 9 s = 1.01 ± 0.02, which should be equal to (1 — (3) with (3 = 0, is recovered. The results for 
the DT model are presented in Fig. 1151 The steady-state persistence exponents have been measured from the power 
law decays shown in Fig. IToT b). In this temporal regime, as shown in panel d), the growth exponent is equal to the 
asymptotic value of 1/5. The persistence behavior of the DT model in this regime is characterized by 9§_ as 0.76 and 
as 0.85, indicating that the relation of Eq. I|14|) holds reasonably well for the (2+1)- dimensional nonlinear MBE 
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FIG. 13: Transient persistence probabilities for the (2+l)-dimensional F and DT growth models. In the case of the F model, 
systems of size L = 1000 x 1000 have been averaged over 200 independent runs, while for the DT model, systems of size 
L = 500 x 500 have been averaged over 800 independent runs. The transient persistence probability for the F model exhibits 
a very fast decay, characterized by a persistence exponent 9 T ~ 6.9 for the last decade of t. More accurate results (see Table 
IVIl are obtained for the exponents in the DT model, although the statistics is not excellent. 



Growth model 


L 


el 


el 


e s + 


e & : 


P 


FM 


200 x 200 


> 6 


> 6 


1.02 ±0.02 


1.00 ±0.02 


0.04 ±0.01 


DT 


200 x 200 


2.84 


2.44 


0.76 ±0.02 


0.85 ±0.02 


0.20 ±0.01 



TABLE V: Transient and steady-state persistence exponents, 9±, for two (2+1) -dimensional discrete growth models. The 
measured value of the growth exponent f3 is shown in the last column. The transient persistence exponents are measured with 
relatively low accuracy due to the rapid temporal decay of the persistence probabilities. 

dynamics, as in the (1 ± l)-dimensional case. It is important to mention that the same values of the persistence 
exponents have been obtained using a DT system with L = 40 x 40, equilibrated for t = 10 5 ML, as required in the 
traditional method (used in most of the (l±l)-dimensional simulations) of measuring the steady-state persistence 
probabilities. Thus the "quick and easy" method of obtaining the steady-state persistence exponent again agrees well 
with the exponent extracted from the actually saturated interface, as discussed above. 

We have also performed some preliminary persistence calculations for the (2±l)-dimensional CKD model in order to 
check the validity of our reported (2±l)-dimensional MBE persistence exponents. Using L — 100 x 100 and to — 1000 
ML, we find that the values of the positive and negative persistence exponents depend to some extent on the chosen 
values for the coefficient A 22 of the nonlinear term and the control parameter c. For example, we obtain 0f « 0.82 and 
9 s . ks 0.77 using A 22 = 5.0 and c = 0.085, and 0f « 0.88 and 0£ w 0.83 using A 22 = 5.0 and c = 0.13. Both cases are 
characterized by a growth exponent of 0.18 ± 0.01, in agreement with Ref. [53j, which is consistent with the expected 
value of 1/5. The results obtained in the latter case are displayed in panel c) of Fig. ^| for the purpose of illustrating 
the similarity between the DT and CKD models. From these observations, we conclude that the (2±l)-dimensional 
DT and CKD persistence results are consistent with each other and they clearly reflect the nonlinearity of the MBE 
dynamical equation in the difference between the values of the positive and negative persistence exponents as expected 
for the up-down asymmetric generic nonlinear situation. 

E. Scaling behavior of the persistence probability 

Since all the results described above have been obtained from simulations of finite systems, it is important to address 
the question of how the persistence probabilities are affected by the finite system size. We have already encountered 
such effects in our study of persistence probabilities for (l±l)-dimensional models (see Table llTlj) . where it was found 
that the measured values of the persistence exponents in the steady state increase slightly as the system size L is 
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FIG. 14: Persistence probabilities for the (2+l)-dimensional F model of size L — 200 x 200, averaged over 200 runs, obtained 
from simulations with different values of the initial time to. a) to — 20 ML. b) to = 200 ML. c) to = 2000 ML. The persistence 
probability curves in case c) show the expected power law decay characterized by the exponent values 9+ = 1.02 ± 0.02 and 
— 1.00 ± 0.02. d) Log-log plot of the interface width W vs deposition time t in units of ML. The slope in the intermediate 
growth decade is f3 ~ 0.04 and thereafter it decreases to zero, as expected. 



decreased, while the value of the growth exponent (3 decreases with decreasing L. We did not investigate finite size 
effects in our study of the transient persistence probabilities because these studies were carried out for large values of 
L and relatively small values of the time t. 

The qualitative dependence of the measured values of the steady-state persistence exponents 6± and the growth 
exponent (3 on the sample size L is not difficult to understand. The steady-state persistence probabilities P±(t 0} h + t) 
exhibit a power law decay with exponent 9± as long as the time t is small compared to the characteristic time scale 
t(L) of the system which is proportional to L z . The decay of P± becomes faster as t approaches and exceeds this 
characteristic time scale. Since this departure from power law behavior occurs at earlier times for smaller systems, 
the value of the persistence exponent extracted from a power law fit to the decay of the persistence probability over 
a fixed time window is expected to increase as the system size is reduced. In a similar way, the measured value of 
(3 is expected to be smaller for relatively small values of L because the precursor to the saturation of the width at 
long times occurs at shorter values of t in smaller systems. Thus, the general trends in the system size dependence 
of the persistence and growth exponents are reasonable. However, it would be useful to obtain a more quantitative 
description of these trends. 

Since the characteristic time scale t(L) ("equilibration" or "saturation" time) of a system of linear size L is 
proportional to L z , one expects, in analogy with the theory of finite size scaling in equilibrium critical phenomena, 
that the steady-state persistence probability P±(t) (in this discussion, we omit the initial time to in the argument of P± 
because the steady-state persistence probability is independent of the choice of t ) would be a function of the scaling 
variable t/L z . Another time scale has to be taken into consideration in a discussion of the scaling behavior of the 
persistence probability. This is the sampling time St which is the time interval between two successive measurements of 
the height at a fixed spatial point. In our simulations of the atomistic growth models, the smallest value of 5t is 1 ML 
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FIG. 15: Persistence probabilities for the (2+l)-dimensional DT model of size L — 200 x 200, averaged over 200 runs, obtained 
from simulations with different values of the initial time io- a) to = 200 ML: The persistence probabilities do no exhibit clear 
power law decay, b) to = 4000 ML: The persistence probability curves show the expected power law decay characterized by 
the exponent values 9+ = 0.77 ± 0.02 and — 0.85 ± 0.02. c) Results for the (2+l)-dimensional CKD model, obtained using 
L = 100 x 100, to — 1000 ML, A = 5 and c = 0.13. d) Log-log plot of the interface width W vs deposition time t in ML for the 
DT model. The slope manifests a crossover from an initial value of ~ 0.26 to the asymptotic value of 0.20. 



because the heights are measured after each deposition of one complete (ideal) monolayer. However, larger integral 
values of St can also be used in the calculation of the persistence probabilities. Since experimental measurements are 
also carried out at discrete time intervals, the presence of a finite value of St has to be accounted for in the analysis 
of experimental data also. Note that the persistence probability itself is mathematically defined, P(to,to + 1), for 
continuous values of time t whereas measurements and simulations are necessarily done on discrete time. 

It has been pointed out in Ref. |5J| that discrete-time sampling of a continuous-time stochastic process does affect 
the measured persistence probability. Such effects have been investigated in detail in the context of a different 
stochastic probability (called the survival probability in Ref. ^(|) that measures the probability of the interface height 
at a fixed position not returning to its time-averaged value within time t. In that work, it was found that the survival 
probability measured for a system of size L with sampling interval 5t is a function of the scaling variables t/L z and 
St/L z . We expect a similar behavior for the steady-state persistence probabilities measured in our simulations. Thus, 
the expected scaling behavior of P±(t,L,5t) is 

Pl(t,L,5t) = f ± (t/L z ,St/L z ), (25) 

_gS 

where the function f±(xi, X2) should decay as x 1 ± for small x\ and X2 <C 1. 

To test the validity of this scaling ansatz, we have carried out calculations of the steady-state persistence probability 
in the linear F model (the positive and negative persistence probabilities are the same in this model) using different 
values of L and St. Due to the linearity of the F model, we have computed a persistence probability P s (t) given by 
the average value of the positive and negative persistence probabilities. If the scaling description of Eq. (|2*51 is valid, 
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FIG. 16: Persistence probability, P(t), for the F model shown for different system sizes with different sampling times. Panel 
a): Double-log plot showing three different persistence curves vs. time corresponding to: L — 100 and St = 4, L — 400 and 
St = 16, L — 800 and St — 81, respectively. Panel b): Finite size scaling of P(t,L,5t). Results for persistence probabilities 
for three different sizes (as in panel a)) with the same value of St/L z (i.e. 1/10 4 ) are plotted versus t/St (z — 2). The dotted 
(dashed) line is a fit of the data to a power law with an exponent of ~ 0.75 (~ 1.0). 



then plots of P s (t, L, St) versus t/L z for different values of L and St should coincide if the value of St for the different 
sample sizes are chosen such that the ratio St/L z remains constant. As shown in Fig. 16 where we present the data 
obtained from simulations of the (l+l)-dimensional F model for three different values (200, 400 and 800) of L and 
three correspondingly different values of St (4, 16 and 64, so that St/L z with z = 2 is held fixed at 10~ 4 ), plots of 
P (t) versus t/5t (which is proportional to t/L z because St is chosen to be proportional to L z ) for the three different 
sample sizes exhibit an excellent scaling collapse. These results confirm the validity of the scaling form of Eq. (|25p. 

As shown in Fig. 16, the scaling function / exhibits the expected power law behavior for relatively small values 
of t/St. Our results also show signatures of a crossover to a power law decay with exponent 1 as t approaches and 
exceeds the characteristic time scale t(L) (this crossover occurs near t/L z ~ 0.1 in the F model). We discuss below 
a possible explanation for this behavior. 

Height fluctuations at times to and t^+t are expected to be completely uncorrelated if t is large compared to t(L). 
Therefore, the persistence probability P (t) for values of t much larger than t(L) may be obtained by considering 
a collection of fluctuating variables which have the same probability distribution (since the system is in the steady 
state), and which are completely uncorrelated with one another. Let xq, xx, x%, . . . represent such a collection of 
variables (these variables may be thought of as the height at a particular site measured at regularly spaced times with 
spacing larger than t(L)). For simplicity, we assume that each Xi is uniformly distributed between —a and a. Then, 
given a particular value of Xq, the probability P + (xq, n) that all the variables x%, 1 < i < n are larger than xq may be 
easily obtained as 



P + (x ,n) 



h 1 ,h 



= [(a-x Q )/(2a)Y 



(26) 



The positive persistence probability P + (n) is obtained by averaging this probability over the probability distribution 
of xq. Thus, we have 



P+(n) 



1 

2a 



P + (x ,n)dx 



1 



n+V 



(27) 



which decays as a power law with exponent 1 for large n. This power law behavior does not depend on the form of 
the assumed probability distribution for the fluctuating variables {x{\. Assuming a general probability distribution 
p(x) with J_ p(x)dx = 1, Eq. (|27|l can be written as 



P+(n) 



p(xo) 



p{x)dx 



dxo 



P(x ) 



p(x)dx 



dxn. 



(28) 



For large n, the quantity that multiplies p(xq) in Eq. (|28|l is of order unity only for values of xq for which J_° p(x)dx 
is of order 1/n. Physically, this means that the positive persistence probability is nonzero for large n only if the initial 
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value xq is very close to the lower limit of the allowed range of values. This effectively restricts the upper limit of the 
integral over xq to yo where yo satisfies the requirement that 



where C is a constant of order unity. Since the quantity that multiplies p(xq) in Eq. (|28|l is of order unity for such 
values of xq, it follows that 



This simple analysis shows that the simulation results for the behavior of the scaling function of Fig. 16 for large 
values of t/St are quite reasonable. 

While we have not carried out similar scaling analyses [Eq. (|25l) ] for other models, we expect the scaling form of 
Eq. (|25(l to be valid in general. We expect that such scaling analysis of the persistence probability as a function of 
the system size L and the sampling time St would be useful in the analysis of numerical and experimental data in 
the future. In fact, we believe that a direct experimental verification of the scaling ansatz defined by Eq. I|25|) will be 
valuable. 



In this paper we have investigated the temporal first passage statistics, expressed in terms of temporal persistence 
probabilities, for a variety of atomistic models that provide discrete realizations of several linear and nonlinear Langevin 
equations for the stochastic dynamics of growing and fluctuating interfaces. Using extensive kinetic Monte Carlo 
simulations, we have obtained transient and steady-state persistence exponents for these (1+1) and (2+l)-dimensional 
SOS and RSOS growth models. We have followed the methodology of Krug et al. |23Ll24| and extended their numerical 
work to the nonlinear MBE dynamical equation by studying the persistence behavior of the atomistic DT, WV, CKD 
and KPK models. From these studies, we have identified two persistence exponents for each of the two temporal 
regimes (transient and steady-state) of the persistence probability. The difference between the values of the two 
exponents reflects the nonlinearity (and the resulting lack of up-down symmetry) of the MBE dynamical equation. 

Among the models studied here, we find that in (1+1) dimensions and in the range of system sizes used in our 
simulations, WV and DT models are hardly distinguishable from the point of view of transient and steady-state 
persistence behavior: the persistence exponents obtained for these two models are very close to each other. We, 
therefore, conclude that in the range of simulation parameters used in this study, the (l+l)-dimensional DT and WV 
models belong to the same universality class (namely the MBE universality class) as far as their persistence behavior is 
concerned. A separate investigation is required in order to understand the universality class of the WV model in (2+1) 
dimensions. The KPK model appears not to reflect well the nonlinear feature of the underlying dynamical equation 
in the values of the positive and negative persistence exponents. This is probably due to strong finite size effects 
arising from the small lattice sizes used in our traditional steady-state simulations (i.e. using to ~ L z ). These finite 
size effects appear to lead to overestimated persistence exponents [and underestimated growth exponent, consistent 
with Eq. ©]. 

We have also investigated the CKD model, which is another discrete model belonging to the MBE universality 
class, our main goal being a closer examination of how the nonlinearity of the underlying continuum equation is 
reflected in values of the transient and steady-state persistence exponents. In this case we have obtained clearly 
different values for the positive and negative persistence exponents. The predictions of the CKD model concerning 
the persistence exponents have been checked by applying the noise reduction technique to the DT model. We found 
that for the MBE universality class, the steady-state persistence exponents in (1+1) dimensions are: 9+ — 0.66 + 0.02 
and = 0.78 ± 0.02. These two values represent the average of the results obtained for the CKD and the noise 
reduced DT models. These results suggest that measurements of persistence exponents can be profitably used to 
detect the presence of nonlinearity in the continuum equations underlying surface growth and fluctuation phenomena. 

The observed difference between the positive and negative steady-state persistence exponents for the models in 
the MBE universality class implies that the relation of Eq. ©, known to be valid for linear Langevin equations 
(our results for the linear F and LC models are in agreement with this relation), can not be satisfied by both these 
exponents. Thus, it is clear that at least one of these steady-state persistence exponents is not related to the usual 
dynamic scaling exponents in a simple way. We have found that the relation of Eq. @ is approximately satisfied 
(within the error bars of our numerical results) by the smaller one of the two steady-state persistence exponents in 
all the (1+1)- and (2+l)-dimensional discrete stochastic growth models studied in this paper. We have also shown 




(29) 




(30) 
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analytically that this relation between the smaller persistence exponent and the growth exponent is, in fact, exact. 
The smaller exponent appears to correspond to the positive (negative) persistence probability if the top (bottom) 
part of the steady-state interface profile is smoother. This observation suggests a deep (and potentially important) 
connection between the surface morphology and the associated persistence exponent, which has no simple analog 
in the dynamic scaling approach where the critical exponents (a, 0, z = a/f3) by themselves do not provide any 
information about the up-down symmetry breaking in the surface morphology. Further investigation of this aspect 
would be very interesting and highly desirable, particularly if experimental information on persistence properties of 
nonequilibrium surface growth kinetics becomes available. 

Our investigation of the effects of the initial configuration on the persistence probabilities indicates that the transient 
persistence exponents can be obtained only if the interface is completely flat at the initial time. This restriction puts 
severe limits on the possibility of measuring the transient persistence exponents in real experiments where it would 
be very difficult, if not impossible, to meet the requirement of zero initial roughness. We have also found the 
surprising and useful result that the steady-state persistence exponents can be accurately measured even if the initial 
configuration corresponds to a value of to that is much smaller than the time required for the interface to reach 
saturation. In other words, the persistence probabilities exhibit their steady-state behavior for measurement times 
comparable to the initial time to even if the value of to is much smaller than L z . This behavior was found in both 
(1+1) and (2+1) dimensions, in all the linear and nonlinear models we studied. This finding is very useful because it 
opens up the possibility of numerically calculating the steady-state persistence exponents for large systems and for 
higher dimensions as well. In fact, this observation enabled us to calculate the steady-state persistence exponents for 
(2+l)-dimensional models belonging to the EW and MBE universality classes. For the MBE universality class, we 
have considered the DT model and found the positive and negative persistence exponents in the steady-state to be 
« 0.76 and « 0.85, respectively, in (2+1) dimensions. 

We have also examined in detail the dependence of the measured steady-state persistence probability in the (1+1)- 
dimensional F model on the sample size L and the sampling interval St which is always finite in simulations and 
experimental measurements. We found that this dependence is described by a simple scaling form. The scaling 
function was found to exhibit power law decay with with exponent 1 for times larger than L z . We have proposed 
a simple explanation for this behavior. We believe that such scaling analysis would prove to be useful in future 
numerical and experimental studies of persistence properties. 

We conclude from the results of this study that persistence probabilities provide a valuable set of tool for investi- 
gating the dynamics of nonequilibrium systems in general, and surface growth and fluctuations in particular. Recent 
experimental studies have shown that the concept of persistence can be applied to analyze the dynamics of fluctu- 
ating steps on Al/Si(lll), Ag(lll) and Pb(lll) surfaces (3^, recorded using STM methods. We believe that in 
view of the importance of thermal and shot-noise fluctuations in the dynamics of growing and fluctuating interfaces, 
theoretical and experimental studies of persistence would play an important role in the analysis of the dynamics of 
nonequilibrium surface growth. 
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